# Replication for study 1, Table 4
# Krawatzek / Soroka
# Defending History? The Impact of Context and Speaker in Russia
# Perspectives on Politics
# DOI

require(nnet)
require(marginaleffects)
require(stargazer)
setwd("Insert your local directory")
load("RegressionDataStudy1.RData")

M1 <- multinom(G2AllGroupsF ~ GR2Group + 
                 GroupCritical + GroupRevisionist + GroupPoliticallyLoyal + GroupNoClearView + C2HistoricalInterest, data=RegressionData)

M2 <- multinom(G2AllGroupsF ~ GR2Group + 
                 GroupCritical + GroupRevisionist + GroupPoliticallyLoyal + GroupNoClearView + C2HistoricalInterest +
                 A1Age + c(A1Age*A1Age) + A2GenderMale + A3Capital + A12SectorEconomyPublicPrivate + A18Education + AWealth , data=RegressionData)

Models <- list(M1, M2)
stargazer(Models, p.auto=FALSE,
          star.cutoffs = c(0.05, 0.01, 0.001),
          column.labels   = c("Punishment 'Yes, it is appropriate'", "Punishment 'No, even more severe'", "No punishment at all"),model.names=TRUE,
          add.lines = list(c("Reference Category", "Don't know", "Don't know", "Don't know")),
          digits=3, type="html", title = "Russia Vignette Experiment Type of Violation",
          out="Table 4.doc")

mfx <- marginaleffects(M1, type = "probs", variables = c("GR2Group", "GroupCritical", "GroupRevisionist", "GroupPoliticallyLoyal", "GroupNoClearView"))
mfxsummary <- summary(mfx)
mfxsummary$estimate <- round(mfxsummary$estimate, digits = 4)
mfxsummary$std.error <- round(mfxsummary$std.error, digits = 4)
mfxsummary$statistic <- round(mfxsummary$statistic, digits = 4)
mfxsummary$p.value <- round(mfxsummary$p.value, digits = 4)
mfxsummary$conf.low <- round(mfxsummary$conf.low, digits = 4)
mfxsummary$conf.high <- round(mfxsummary$conf.high, digits = 4)
write.xlsx(mfxsummary, "AME Study 1 without controls.xlsx")

mfx <- marginaleffects(M2, type = "probs", variables = c("GR2Group", "GroupCritical", "GroupRevisionist", "GroupPoliticallyLoyal", "GroupNoClearView"))
mfxsummary <- summary(mfx)
mfxsummary$estimate <- round(mfxsummary$estimate, digits = 4)
mfxsummary$std.error <- round(mfxsummary$std.error, digits = 4)
mfxsummary$statistic <- round(mfxsummary$statistic, digits = 4)
mfxsummary$p.value <- round(mfxsummary$p.value, digits = 4)
mfxsummary$conf.low <- round(mfxsummary$conf.low, digits = 4)
mfxsummary$conf.high <- round(mfxsummary$conf.high, digits = 4)
write.xlsx(mfxsummary, "AME Study 1 with controls.xlsx")